
# ------------------------------------------------------------------------------
# Replication Materials
# 
# title: Coalition Size, Direct Democracy, and Public Spending
# journal: Journal of Public Policy
# authors: Patrick Emmenegger, Lucas Leemann, and André Walter
# date: Sept 2020
# ------------------------------------------------------------------------------




# Main Analysis









m1 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0+I(finex==0)+index_law_ref+numpar+sec_sh+
           fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix+fisrul, findat_oLG)
cov1 <- coeftest(m1, cluster.vcov(m1, cluster = ~canton))[,2]
pval1 <- coeftest(m1, cluster.vcov(m1, cluster = ~canton))[,4]


m2 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init+refmon_adj_imp0+I(finex==0)+index_law_ref*numpar+numpar+sec_sh+
           fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix+fisrul, findat_oLG)
cov2 <- coeftest(m2, cluster.vcov(m2, cluster = ~canton))[,2]
pval2 <- coeftest(m2, cluster.vcov(m2, cluster = ~canton))[,4]

m3 <- lm(dlnexp_pc~factor(canton)+factor(year)+llnexp_pc+index_law_init+numpar*refmon_adj_imp0+I(finex==0)+index_law_ref+sec_sh+
           fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix+fisrul, findat_oLG)
cov3 <- coeftest(m3, cluster.vcov(m3, cluster = ~canton))[,2]
pval3 <- coeftest(m3, cluster.vcov(m3, cluster = ~canton))[,4]


m4 <- lm(dlnexp_pc~factor(canton)+factor(year)+numpar+llnexp_pc+index_law_init*numpar+refmon_adj_imp0*numpar+I(finex==0)+index_law_ref*numpar+sec_sh+
           fir_sh+deprat+babydeathratio_impute+left+PR+phydens+lnpop+lnfed_mix+fisrul, findat_oLG)
cov4 <- coeftest(m4, cluster.vcov(m4, cluster = ~canton))[,2]
pval4 <- coeftest(m4, cluster.vcov(m4, cluster = ~canton))[,4]



# current table 1
texreg(list(m1,m2,m3,m4), omit.coef = "(canton|year)", booktabs = T, dcolumn = T, 
       stars = c(0.001, 0.01, 0.05, 0.1), override.se = list(cov1,cov2,cov3,cov4),
       override.pvalues = list(pval1,pval2,pval3,pval4),
       use.packages = F, scalebox = .8, label = "est1",caption = "Direct Democracy and Government Spending, 1860-2015",
       file = "out/est1.tex", caption.above = TRUE,float.pos="h!",
       custom.coef.map = list("numpar" = "Number Parties","index_law_init" = "Initiatives",
                              "refmon_adj_imp0" = "Financial Ref. (Threshold)",
                              "I(finex == 0)TRUE" = "With Financial Referednum",
                              "index_law_ref" = "Law Ref.","llnexp_pc" = "Lag Dep. Var.",
                              "fir_sh" = "Share First Sector","sec_sh" = "Share Second Sector",
                              "deprat" = "Dependency Ration", "babydeathratio_impute" = "Infant Mort.",
                              "left" = "Share Left Parties","PR" = "Proportional Representation",
                              "phydens" = "Physician Density",
                              "lnpop" = "ln Population Size","lnfed_mix" = "ln Federal Subsidies",
                              "numpar:index_law_init" = "Num. Par. * Initiative",
                              "numpar:index_law_ref" = "Num. Par. * Law Referednum",
                              "numpar:refmon_adj_imp0" = "Num. Par. * Fin. Referendum",
                              "numpar:finex" = "Num. Par. * Existence of Fin. Referendum"))


###############################################################################################################
###############################################################################################################




######################################################################
################################# Plots ##############################

#### marginal Effect Initiative
N.sim <- 1000
set.seed(123)
BETA <- mvrnorm(N.sim,coef(m4), vcov(m4))

BETA.sel <- BETA[,c("index_law_init","numpar:index_law_init")]

ME.init <- matrix(NA,N.sim,5)

ME.init[,1] <- BETA.sel[,1] + 1*BETA.sel[,2]
ME.init[,2] <- BETA.sel[,1] + 2*BETA.sel[,2]
ME.init[,3] <- BETA.sel[,1] + 3*BETA.sel[,2]
ME.init[,4] <- BETA.sel[,1] + 4*BETA.sel[,2]
ME.init[,5] <- BETA.sel[,1] + 5*BETA.sel[,2]

ME.init1 <- matrix(NA,5,5)

for (i in 1:5){
  ME.init1[,i] <- quantile(ME.init[,i], c(0.025,0.05,.5,0.95, 0.975))
}

x.init1 <- matrix(NA,5,5)
for (i in 1:5){
  x.init1[,i] <- i
}



pdf("out/ME_init.pdf",width=9,height=5.5)

par(family="CMU Serif")
plot(x.init1[3,],ME.init1[3,], bty="n", ylim=c(-.15,.1),
     ylab="Marginal Effect of Initiative", xlab="Number of Political Parties",
     pch=19, col="blue", yaxt='n')
axis(2,c(-.1,-.05,0,.05))

for (i in 1:5){
  segments(x.init1[2,i],ME.init1[2,i],x.init1[4,i],ME.init1[4,i], lwd=4, col="blue")  
  segments(x.init1[2,i],ME.init1[1,i],x.init1[4,i],ME.init1[5,i], lwd=1, col="blue")  
}  
abline(h=0, lty=2)

c.part <- table(findat$numpar)/(10*sum(table(findat$numpar)))
for (i in 1:5){
  rect(i-0.1,-.15,i+0.1,-.15+c.part[i],col=rgb(0,0,255,100,maxColorValue = 255), border = NA)
}



mtext("Share of \n Cases", side = 2, adj = .01, col= "blue")

dev.off()



# #### marginal Effect Referendum
N.sim <- 1000
set.seed(123)
BETA <- mvrnorm(N.sim,coef(m4), vcov(m4))

BETA.sel <- BETA[,c("index_law_ref","numpar:index_law_ref")]

ME.init <- matrix(NA,N.sim,5)

ME.init[,1] <- BETA.sel[,1] + 1*BETA.sel[,2]
ME.init[,2] <- BETA.sel[,1] + 2*BETA.sel[,2]
ME.init[,3] <- BETA.sel[,1] + 3*BETA.sel[,2]
ME.init[,4] <- BETA.sel[,1] + 4*BETA.sel[,2]
ME.init[,5] <- BETA.sel[,1] + 5*BETA.sel[,2]

ME.init1 <- matrix(NA,5,5)

for (i in 1:5){
  ME.init1[,i] <- quantile(ME.init[,i], c(0.025,0.05,.5,0.95, 0.975))
}

x.init1 <- matrix(NA,5,5)
for (i in 1:5){
  x.init1[,i] <- i
}



pdf("out/ME_ref.pdf",width=9,height=5.5)

par(family="CMU Serif")
plot(x.init1[3,],ME.init1[3,], bty="n", ylim=c(-.15,.09),
     ylab="Marginal Effect of Law Referendum", xlab="Number of Political Parties",
     pch=19, col="blue", yaxt='n')
axis(2,c(-.1,-.05,0,.05))

for (i in 1:5){
  segments(x.init1[2,i],ME.init1[2,i],x.init1[4,i],ME.init1[4,i], lwd=4, col="blue")
  segments(x.init1[2,i],ME.init1[1,i],x.init1[4,i],ME.init1[5,i], lwd=1, col="blue")
}
abline(h=0, lty=2)

c.part <- table(findat$numpar)/(10*sum(table(findat$numpar)))
for (i in 1:5){
  rect(i-0.1,-.15,i+0.1,-.15+c.part[i],col=rgb(0,0,255,100,maxColorValue = 255), border = NA)
}

mtext("Share of \n Cases", side = 2, adj = .01, col= "blue")

dev.off()

# #### marginal Effect Financial Referendum
N.sim <- 1000
set.seed(123)
BETA <- mvrnorm(N.sim,coef(m4), vcov(m4))

#BETA.sel <- BETA[,c("refmon_adj4","numpar:refmon_adj4")]
BETA.sel <- BETA[,c("refmon_adj_imp0","numpar:refmon_adj_imp0")]

ME.init <- matrix(NA,N.sim,5)

ME.init[,1] <- BETA.sel[,1] + 1*BETA.sel[,2]
ME.init[,2] <- BETA.sel[,1] + 2*BETA.sel[,2]
ME.init[,3] <- BETA.sel[,1] + 3*BETA.sel[,2]
ME.init[,4] <- BETA.sel[,1] + 4*BETA.sel[,2]
ME.init[,5] <- BETA.sel[,1] + 5*BETA.sel[,2]

ME.init1 <- matrix(NA,5,5)

for (i in 1:5){
  ME.init1[,i] <- quantile(ME.init[,i], c(0.025,0.05,.5,0.95, 0.975))
}

x.init1 <- matrix(NA,5,5)
for (i in 1:5){
  x.init1[,i] <- i
}

pdf("out/ME_Fref.pdf",width=9,height=5.5)

par(family="CMU Serif")
plot(x.init1[3,],ME.init1[3,], bty="n", ylim=c(-.8,0.6),
     ylab="Marginal Effect of Financial Referendum", xlab="Number of Political Parties",
     pch=19, col="blue", yaxt='n')
axis(2,c(-0.2,0,.2,.4,.6,.8,1.0,1.2))

for (i in 1:5){
  segments(x.init1[2,i],ME.init1[2,i],x.init1[4,i],ME.init1[4,i], lwd=4, col="blue")
  segments(x.init1[2,i],ME.init1[1,i],x.init1[4,i],ME.init1[5,i], lwd=1, col="blue")
}
abline(h=0, lty=2)

c.part <- table(findat$numpar)/(2*sum(table(findat$numpar)))
for (i in 1:5){
  rect(i-0.1,-.85,i+0.1,-.85+c.part[i],col=rgb(0,0,255,100,maxColorValue = 255), border = NA)
}

mtext("Share of \n Cases", side = 2, adj = .01, col= "blue")

dev.off()





###############################################################################################################
###############################################################################################################

#### Plot Longterm Effects


N.sim <- 10000
set.seed(987)
BETA <- mvrnorm(N.sim,coef(m4), vcov(m4))



pdf("out/LT_init1.pdf",width=17,height=7)

init <- 1.5
partyL <- 1
a <- (exp((BETA[,"index_law_init"]*init + BETA[,"numpar:index_law_init"]*partyL*init)/(-BETA[,"llnexp_pc"]))-1)*100

partyL <- 5
b <- (exp((BETA[,"index_law_init"]*init + BETA[,"numpar:index_law_init"]*partyL*init)/(-BETA[,"llnexp_pc"]))-1)*100
QQ1 <- quantile(a-b, c(0.05,.5,.95))



init <- 5
partyL <- 1
a2 <- (exp((BETA[,"index_law_init"]*init + BETA[,"numpar:index_law_init"]*partyL*init)/(-BETA[,"llnexp_pc"]))-1)*100

partyL <- 5
b2 <- (exp((BETA[,"index_law_init"]*init + BETA[,"numpar:index_law_init"]*partyL*init)/(-BETA[,"llnexp_pc"]))-1)*100
QQ2 <- quantile(a2-b2, c(0.05,.5,.95))



amou1 <- c(min(a-b) - 0.15*diff(range(a-b)),max(a-b) + 0.15*diff(range(a-b)))
amou2 <- c(min(a2-b2) - 0.15*diff(range(a2-b2)),max(a2-b2) + 0.15*diff(range(a2-b2)))
par(mfrow=c(1,2), family="CMU Serif",mar=c(5,5,2,0))
hist(a-b, breaks=seq(floor(amou1[1]),ceiling(amou1[2]),length.out = 100), col=rgb(0,0,255,100, maxColorValue = 255), border="white",
     xlab=" ", ylab="Low Levels of Initiative Rights", xlim=amou1,  cex.lab=2, main="")
abline(v=QQ1[1], col="blue", lty=2)
abline(v=QQ1[3], col="blue", lty=2)
abline(v=QQ1[2], col="blue", lty=1, lwd=3)

hist(a2-b2, breaks=seq(floor(amou2[1]),ceiling(amou2[2]),length.out = 100), col=rgb(255,0,0,100,maxColorValue = 255), border="white",
     xlab="", ylab="High Levels of Initiative Rights", xlim=amou2, cex.lab=2, main="")
abline(v=QQ2[1], col="red", lty=2)
abline(v=QQ2[3], col="red", lty=2)
abline(v=QQ2[2], col="red", lty=1, lwd=3)

dev.off()

sink("out/LR_1.tex")
cat(paste(as.character(round(QQ1[2], digits=1)),"%", sep=""), "\n")
sink()


sink("out/LR_2.tex")
cat(paste(as.character(round(QQ2[2], digits=1)),"%", sep=""), "\n")
sink()










pdf("out/LT_ref1.pdf",width=17,height=7)
init <- 1
partyL <- 2
a <- (exp((BETA[,"index_law_ref"]*init + BETA[,"numpar:index_law_ref"]*partyL*init)/(-BETA[,"llnexp_pc"]))-1)*100

partyL <- 4
b <- (exp((BETA[,"index_law_ref"]*init + BETA[,"numpar:index_law_ref"]*partyL*init)/(-BETA[,"llnexp_pc"]))-1)*100
QQ1 <- quantile(b-a, c(0.05,.5,.95))



init <- 6
partyL <- 2
a2 <- (exp((BETA[,"index_law_ref"]*init + BETA[,"numpar:index_law_ref"]*partyL*init)/(-BETA[,"llnexp_pc"]))-1)*100

partyL <- 4
b2 <- (exp((BETA[,"index_law_ref"]*init + BETA[,"numpar:index_law_ref"]*partyL*init)/(-BETA[,"llnexp_pc"]))-1)*100
QQ2 <- quantile(b2-a2, c(0.05,.5,.95))



amou1 <- c(min(b-a) - 0.15*diff(range(b-a)),max(b-a) + 0.15*diff(range(b-a)))
amou2 <- c(min(b2-a2) - 0.15*diff(range(b2-a2)),max(b2-a2) + 0.15*diff(range(b2-a2)))
par(mfrow=c(1,2), family="CMU Serif",mar=c(5,5,2,0))
hist(b-a, breaks=seq(floor(amou1[1]),ceiling(amou1[2]),length.out = 100), col=rgb(0,0,255,100, maxColorValue = 255), border="white",
     xlab=" ", ylab="Low Levels of Law Referendum Rights", xlim=amou1,  cex.lab=2, main="")
abline(v=QQ1[1], col="blue", lty=2)
abline(v=QQ1[3], col="blue", lty=2)
abline(v=QQ1[2], col="blue", lty=1, lwd=3)

hist(b2-a2, breaks=seq(floor(amou2[1]),ceiling(amou2[2]),length.out = 100), col=rgb(255,0,0,100,maxColorValue = 255), border="white",
     xlab="", ylab="High Levels of Law Referendum Rights", xlim=amou2, cex.lab=2, main="")
abline(v=QQ2[1], col="red", lty=2)
abline(v=QQ2[3], col="red", lty=2)
abline(v=QQ2[2], col="red", lty=1, lwd=3)

dev.off()




sink("out/LR_3.tex")
cat(paste(as.character(round(QQ1[2], digits=1)),"%", sep=""), "\n")
sink()


sink("out/LR_4.tex")
cat(paste(as.character(round(QQ2[2], digits=1)),"%", sep=""), "\n")
sink()


# Fref Plot
N.sim <- 1000
set.seed(123)
BETA <- mvrnorm(N.sim,coef(m1) , vcov(m1) )

pdf("out/LT_Fref1.pdf",width=17,height=7)

#Fref <- quantile(findat2$refmon_adj/100, c(0.25,0.75), na.rm = TRUE)
Fref <- quantile(findat2$refmon_adj[findat2$year==1930]/100, c(0.25,0.75),na.rm=TRUE)

a <- (exp(BETA[,"refmon_adj_imp0"]*Fref[1]/(-BETA[,"llnexp_pc"]))-1)*100
QQ <- quantile(a, c(0.025,.5,.975))
QQ

b <- (exp(BETA[,"refmon_adj_imp0"]*Fref[2]/(-BETA[,"llnexp_pc"]))-1)*100
QQ1 <- quantile(b, c(0.025,.5,.975))
QQ1

QQ2 <- quantile(a-b, c(0.025,.5,.975))


amou1 <- c(-1,7)
amou2 <- c(-2,13)
amou3 <- c(-8,2)


par(mfrow=c(1,3), family="CMU Serif",mar=c(5,5,2,0))
hist(a, breaks=seq(floor(min(a)),ceiling(max(a)),length.out = 50), col=rgb(0,0,255,100,maxColorValue = 255), border="white",
     xlab=" ", ylab="High Threshold", xlim=amou1, main="", cex.lab=2, cex.main=2)
abline(v=QQ[1], col="blue", lty=2)
abline(v=QQ[3], col="blue", lty=2)
abline(v=QQ[2], col="blue", lty=1, lwd=3)


hist(b, breaks=seq(floor(min(b)),ceiling(max(b)),length.out = 50), main="", col=rgb(255,0,0,100,maxColorValue = 255), border="white",
     xlab="", ylab="Low Threshold", xlim=amou2, cex.lab=2, cex.main=2)
abline(v=QQ1[1], col="red", lty=2)
abline(v=QQ1[3], col="red", lty=2)
abline(v=QQ1[2], col="red", lty=1, lwd=3)


hist(a-b, breaks=seq(floor(min(a-b)),ceiling(max(a-b)),length.out = 50), main="", col=rgb(255,0,255,100,maxColorValue = 255), border="white",
     xlab="", ylab="Difference in Long-Run Effect", xlim=amou3, cex.lab=2, cex.main=2)
abline(v=QQ2[1], col="purple", lty=2)
abline(v=QQ2[3], col="purple", lty=2)
abline(v=QQ2[2], col="purple", lty=1, lwd=3)


dev.off()



sink("out/LR_5.tex")
cat(paste(as.character(round(QQ2[2], digits=1)),"%", sep=""), "\n")
sink()